On the ground— state energy of finite Fermi systems 
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We study the ground-state shell correction energy of a fermionic gas in a mean-field approxima- 
tion. Considering the particular case of 3D harmonic trapping potentials, we show the rich variety of 
different behaviors (erratic, regular, supershells) that appear when the number-theoretic properties 
of the frequency ratios are varied. For self-bound systems, where the shape of the trapping potential 
is determined by energy minimization, we obtain accurate analytic formulas for the deformation and 
the shell correction energy as a function of the particle number TV. Special attention is devoted to 
the average of the shell correction energy. We explain why in self-bound systems it is a decreasing 
(and negative) function of TV". 

PACS numbers: 21.10.Dr, 03.65.Sq, 21.60.-n 



I. INTRODUCTION 



When estimating the ground-state energy of interact- 
ing Fermi gases, one can schematically rely on two dif- 
ferent approaches. For few-particle systems (say, less 
than 10-20), direct ab initio calculations are available 
at present For larger systems, and in particular for 
global mass calculations in atomic nuclei, approximate 
schemes based either on mean-field calculations [2[ or 
shell models 0, H| are inevitable. In mean field theories, 
in which we focus here, the energy or mass of the sys- 
tem is naturally decomposed into a smooth part and a 
fluctuating part. Methods like Thomas-Fermi (TF) or 
Wigner-Kirkwood theories provide an approximation to 
the smooth part whereas the oscillatory shell struc- 

ture may be described by semiclassical methods [4]. In 
practice, these two contributions may be well separated 
in a grand canonical scheme, where the energy is consid- 
ered as a function of the chemical potential \x: the TF 
contribution describes the smooth dependence on chemi- 
cal potential, whereas the oscillatory component has zero 
average (with respect to fi) and describes deviations with 
respect to the mean behavior. 

In isolated systems with a well defined number of par- 
ticles TV, things are different. When the dependence of 
the ground-state energy E(N) is considered as a function 
of TV, one can show [5[ that its fluctuating part, E(N), 
in contrast to usual expectations, has a non-zero aver- 
age (as a function of TV). It follows that the fluctuating 
part contributes to the smooth part E(N) of the energy, 
and the frontier between smooth and fluctuating com- 
ponents is blurred. This is a generic effect, although its 
importance depends on the symmetries (and intrinsic dy- 
namics) of the mean field potential. A description of this 
effect was recently developed [j| for a confining poten- 
tial which keeps its shape fixed (up to a possible scaling 
factor) when the number of particles varies (cf Eq. ([6|) be- 



low). 

Although the latter situation may be relevant in many 
experimental set-ups, like cold dilute Fermi gases in mag- 
netic atom traps [6|, where the external HO potential 
dominates over the mean-field interaction energy Q , an- 
other relevant case is that of self-bound systems, like the 
atomic nucleus. In these systems, an effect that was not 
included in the previous description appears: the shape 
of the average self-bound confining potential depends on 
the number of particles. As is well known, at a given TV 
the shape is determined by minimizing the energy of the 
system. The minimization of the smooth part of the en- 
ergy leads generically to an isotropic shape. We denote 
this contribution E sp h(N). For that shape, and for par- 
ticular values of TV (magic numbers) the contribution of 
the fluctuating part E sp h(N) is large and negative, thus 
reducing the total energy with respect to E sp h (N) . Away 
from magic numbers, the amplitude of E sp h(N) rapidly 
decreases, and may eventually become positive. In order 
to avoid this behavior, the system deforms, trying to keep 
the value of the oscillating part of the self-bound energy 
E S b(N) negative and as large as possible. Though in 
most realistic situations E SB (N) <C E(N), the behavior 
of E SB (N) has a strong influence on the shape of the sys- 
tem. Schematically (see Fig|4j, the behavior of E SB (N) 
is therefore a fluctuating negative function of TV, with 
larger amplitudes around the magic numbers. 

It follows that in self-bound systems, as for Fermi 
gases confined by a fixed external potential, the aver- 
age part of the energy fluctuations is again generically 
different from zero. However, the properties of the av- 
erage are very different in these two cases. In partic- 
ular, the average is positive for a fixed shape, while it 
is negative in self-bound systems. Such a bias toward 
negative energies of the fluctuating part in self bound 
systems is clearly observed in realistic calculations 
In the bottom part of Fig. [5] we plot the nuclear ground- 
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state shell correction energy computed in Ref.[8| using 
a macroscopic-microscopic model, whose results are in 
good agreement with experimental data. We do observe 
a tendency toward negative values, with a non-zero slope 
for the average part of the fluctuations. One of our pur- 
poses here is to provide a quantitative description of this 
effect. 

We will consider the case of non-interacting fermions 
in 3D whose self-bound confining potential is assumed to 
have a quadratic harmonic shape. The reason for such a 
choice is purely technical, since it allows, to some extent, 
for an explicit analytic description. In spite of its sim- 
plicity, we will show that it provides a correct description 
of what is observed in more realistic calculations. 

Harmonic potentials were intensively investigated in 
the past, starting from the Nilsson model of the nuclear 
deformation [§] . This is an integrable (separable) model. 
It is known that the statistical properties of its single- 
particle spectrum do not coincide with the generic (Pois- 
son) behavior expected in integrable systems [lcj . The 
statistics, explored mainly in D=2, are in fact very sensi- 
tive to the number theoretic properties of the frequency 
ratios we will show these number theo- 

retic properties also strongly influence the behavior of 
the many body system. 

Our analysis of the minimization of the energy of the 
Fermi gas will be based on a semiclassical theory. Though 
this theory is exact for the single-particle density of 
states, the different cases (of irrational or rational fre- 
quency ratios) should, however, be considered carefully 
ll| . We shall see that the amplitude and phase of shell 
effects are directly controlled by the number-theoretic 
properties of the frequencies. The output of the min- 
imization problem thus depends on a delicate interplay 
between these number-theoretic properties and the num- 
ber of particles in the gas. 

The manuscript is organized as^follows: in section II 
we study the main properties of E(N) for an harmonic 
potential of given frequency ratios. We will illustrate 
the large variety of different behaviors that could emerge 
(including shell and supershell structures). In section III, 
we deal with the case of self-bound systems, and apply 
our results to a schematic model of the atomic nucleus 
(that we compare to realistic calculations). 



I. GENERAL SETTING 

For a given potential the single-particle level density 
is decomposed into a smooth part p coming from the TF 
theory plus an oscillatory contribution p: 



p(e) = p(e) + p(e) . 



(1) 



Here, e is the single-particle energy. Below, we will show 
in more detail how to describe these two components. 



The spin degeneracy is included in the level density. 
Equation |T]) induces a corresponding decomposition of 
the integrated level density (that counts the number of 
single-particle states up to an energy p): 



M(p) =Af(p)+M(p) 



p(e) de 



p{e)de. (2) 



In order to study a system with a finite number of par- 
ticles, we use canonical expressions for thermodynamic 
quantities. For a system of TV non-interacting fermions, 
we define its ground state energy E(N), the shell correc- 
tion energy E(N) and the smooth TF component E(N) 
as[|i: ' 

E(N) = E(N) - E(N) = [ " ep{e)de - [ " ep{e)de. 

Jo Jo 

(3) 

The chemical potential p N and its smooth part p N fix the 
number of particles. They are defined by inversion of the 
exact and average integrated level densities, Af(p N ) = N 
and 



(4) 



respectively. Unfortunately, Eq.([3]) is difficult to exploit 
analytically because the discretization of p N is difficult 
to impose. From Eq.Q it can be shown that, neglecting 
terms of second order in the parameter p — p, E may be 
approximated by 12|, [I3J : 



E{N) 



M(e)de 



(5) 



This, together with the definition of E(N), are the ba- 
sic equations upon which our analysis of ground state 
energies of Fermi gases will be based on. 

If, in Eq. ([5]), p N is considered as a continuous vari- 
able, it can easily be shown that it gives a wrong result 
for the average of the fluctuating function E(N). In- 
deed in a system with a fixed number of particles, the 
chemical potential takes discrete values as TV varies. The 
fluctuating part of the ground state energy is sampled 
at particular values of the chemical potential, implying 
a modification of its average value. Recently, an explicit 
description for this effect was given. It was found that 
the contribution of the fluctuating part to the average of 
the energy is given by || 



(E(N)) N ^(Af 2 (p N ))n /p{p h 



(6) 



where we use brackets to denote an average over an ap- 
propriate chemical potential window and to distinguish 
this contribution with respect to the TF smooth term. 
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II. EXTERNAL POTENTIALS 
A. Triaxial case with irrational frequencies 

We consider a particle in a 3D harmonic potential. The 
Hamiltonian is given by: 

H = \tpl+P 2 y +pl) + \(Qlx 2 + Qly 2 + Qlz 2 ) , (7) 

where all quantities are dimensionless (units h = m = 
1). Qi, Qi and Q3 are the frequencies of the harmonic 
oscillator (HO). In this subsection we considered them as 
incommensurable real numbers with pairwise irrational 
ratios. Using Eq.Q and the expression of the TF level 
density p given in Table I, we calculate from Eq.|4|) the 
smooth chemical potential to leading order in N: 

= (3QiQ 2 Q 3 A r ) 1/3 • (8) 

For irrational frequency ratios, the only periodic orbits 
of the system are the one-dimensional oscillations along 
the three principal axis x, y and z. These orbits are iso- 
lated, and are the backbone for the description of the 
fluctuating part of the different quantities. In particular, 
the single particle level density p(e) was computed in (llj 
(using semiclassical methods based on the trace formula, 
which coincide with exact results obtained from the in- 
verse Laplace transform of the exact partition function) . 
From Eqs.([2]) and $5§, and the previous expression ([5]) 
of p N , we compute analytically the expressions for A/"(/i) 
and E(N). These, together with that of p (see Ref . (Til] ) . 
are reported in Table II (first column, case A). 

From periodic orbit theory (POT), it is known that 
the main features of the shell correction energy can be 
recovered by taking into account only the first (short- 
est) orbits [13| . which correspond to the first terms in 
each of the three sums (lowest k values) of E(N). As 
the number of particles N is varied, the interferences be- 
tween these terms produce an oscillatory pattern. For 
potentials in which the single-particle spectrum has a 
simple structure, like the isotropic HO discussed in the 
next subsection, the oscillatory pattern of the energy as 
a function of N will be quite regular. One speaks in that 
case of shell effects and magic numbers (which correspond 
to the minima of the regular oscillation). In the present 
case of irrational frequency ratios, the structure of the 
single-particle spectrum is highly non-trivial, and leads 
to a complicated pattern of oscillations of the energy as 
a function of N. This is illustrated for a particular case 
in Fig. [TJ The full line is obtained from the correspond- 
ing equation of Table II, using only the first term in each 
sum (k — 1). Although details are missing, note the good 
overall agreement obtained with just three orbits. Note 
the absence of regularity of the pattern. 

As the classical orbits for irrational frequencies do not 
form families but are instead isolated, the amplitude of 



the fluctuations is small (compared to the results of, e.g., 
the isotropic case of the next subsection, see Fig. [2]). The 
purely quantum mechanical counterpart of this state- 
ment is that due to the absence of symmetries in the 
system, no systematic degeneracies occur in the single- 
particle spectrum, fluctuations are small, and accordingly 
the energy does not deviate significantly from its average 
part. Moreover, the typical value of Af is also small, and 
we find a correction to the mean value coming from the 
oscillatory part, calculated from Eq.©, that vanishes as 
jy-2/3 |gj f or f ur ther details of the method). 




FIG. 1: (Color online) Ground state shell correction energy as 
a function of N for Qi — v2, Q2 = v3 and Q3 = \/5. Quan- 
tum computation in dots, theoretical prediction (see Table II) 
using the first term in each sum in full line. 



B. Isotropic case 

Let's consider now the case Qi = Q2 = Qs = L The 
system possesses the U(3) symmetry. Within POT, a 
perturbative approach of the irrational case allows to 
compute the level density [ll|, to obtain the well known 
result given in the top row of the second column of Ta- 
ble II. To compute Af(p,) and E(N), we have followed 
the method explained above (see the results in the sec- 
ond column of Table II). The high degree of symmetry 
of the system leads in this case to a single characteristic 
period for the periodic orbits, implying a much more reg- 
ular pattern of the shell oscillations (compare the upper 
part of Fig. [2] with Fig. Q]). With respect to the triax- 
ial irrational case, the amplitude of the fluctuations is 
much larger (this is related to the high degeneracy of the 
energy levels or, semiclassically, to the fact that the peri- 
odic orbits form families). The approximate frequency of 
the shell fluctuations is given by the phase of the cosine 
function of the first k = 1 term of the sum in E(N). Thus 
magic numbers, given by the values of N that minimize 
E(N), are well approximated by 

i 3 

N MAGIC = T , where i e N* , (9) 
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FIG. 2: (Color online) Upper panel: Ground state shell cor- 
rection energy as a function of N for the isotropic HO. Quan- 
tum computation in dashed line, semiclassical theoretical pre- 
diction in full line. The dotted line correspond to the analyt- 
ical expression of the minima of E (see last row of Table II) . 
Lower panel: Numerical average of the exact fluctuating func- 
tion of the upper panel (dashed line), compared to the theo- 
retical prediction, Eq.© (full line) (from Ref.[5j]). 



which corresponds to integer values of the chemical po- 
tential ft N . 

Classically, all the orbits are closed (periodic), and 
have the same period. This high degeneracy of the peri- 
odic orbits is reflected in the prefactor TV 2 / 3 in front of 
E(N) (cf Table II), as compared to the case of irrational 
ratios. Because there is only one characteristic frequency 
associated to the periodic orbits, no beating effects are 
observed, and the fluctuations are very regular (form- 
ing "shells"). In particular, no supershell structure is 
present, like the characteristic beating pattern observed 
in a spherical cavity with hard walls, produced by the in- 
terference between the triangular and the square orbits, 
see Refe.llEJ. 

In Fig. [51 we clearly see that the average of E(N) is non 
zero and positive, and increases as N increases. Using in 
Eq.([6]) the corresponding expression for J\f(p, N ), we get an 
analytical expression for (E(N)) N (cf fourth row of Table 
II, and [5j for further details). Theory and numerics are 
compared in the bottom part of that figure. 



C. Triaxial case with rational frequencies 

Now Qi, Q2 and Q3 are positive integers and have 
pairwise irreducible ratios. 

Since frequencies are integer, it produces families of 
classical periodic orbits (Lissajous figures). A pertur- 
bative treatment, similar to the isotropic case, can be 
done to obtain an explicit expression for the level den- 
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FIG. 3: (Color online) Upper panel: Ground state shell cor- 
rection energy as the function of N 1 ^ 3 for Q\ = 17, Q2 = 18 
and Q3 = 19. In dots, numerical computation. We observe 
triaxial to isotropic transition for 1Z T ^ 1, which corresponds 
to iV c = 131060 (or 7V C 1/3 - 51). 

Lower left panel: Magnification of the upper panel for low 
values of iV (1Z T > 7). In dots, numerical computation. 
The full line corresponds to the analytical expression of E T 
with the first term in each sum. The triaxial component of E 
dominates, and shows super shell effects due to interferences 
of the shortest orbits. 

Lower right panel: Magnification of the upper panel for 
iV 1/3 - 150 {TZt ~ 0.11). Dots, numerical computa- 
tion. The full line corresponds to the analytical expression of 
Ei. The isotropic component of E dominates. The behavior 
is similar to the isotropic case of Fig. [2] 



sity [llj. For convenience, the oscillating part of the level 
density is decomposed into two different parts (see Table 
II): f>j (resp. p T ) which is connected to the level density 
of the isotropic (resp. triaxial with irrational frequency 
ratios) HO. Below, these two components are referred to 
as " isotropic" and " triaxial" , respectively. We follow the 
same method as before in order to get expressions for 
Kf{jj) and E(N) (cf third column of Table II). 

To understand in this case the TV-dependence of the 
fluctuating part of the energy, it is useful to calculate the 
ratio 



^ At 



(10) 



of the typical amplitude of the triaxial component with 
respect to the isotropic one. These two amplitudes are 
approximated here by computing the amplitude of the 
first term k = 1 of the corresponding sum (cf Table II) 



A, 



(3A0 



2/3 



27r 2 (Q 1 Q 2 Q 3 ) 1 /3 



/ Qi/(47r 2 |sm(7rQ2/Qi)sin(7rQ 3 /Qi) 
max Q 2 /(47r 2 | sin(7rQi/Q 2 ) sin(7r(33/(32) 
V Q3/(47r 2 |sin(7rQi/Q 3 )sin(7rQ2/Q3) 



(11) 
(12) 
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A finer estimate of A T requires an analysis of the num- 
ber theoretic properties of the k-th dependence of the 
denominators in E T , that we shall not do here. Since, 
for given frequencies, A T is a constant independent of 
N, 1Z T oc A~ 2 / 3 . Thus, for any rational set of fre- 
quencies, there exists a critical number of particles N c 
above which the behavior of E(N) is dominated by the 
isotropic component. In this high- A regime the behav- 
ior of E(N) qualitatively coincides, up to a rescaling of 
the overall amplitude by (QiQ2Qs)~ 1 ^ 3 and of a rescal- 
ing of the phase by (<9iQ2<33) 1//3 , with the behavior of 
the isotropic HO described in the subsection II. B. This 
corresponds to a regular (i.e. single-frequency) pattern 
of large amplitude oscillations (or shells), with clearly 
defined magic numbers, as a function of N (compare 
Fig. [5] and the lower right panel of Fig. [3]). The value 
of the magic numbers may be estimated similarly to the 
previous subsection. This behavior illustrates the occur- 
rence of an increasing number of accidental degeneracies 
in the single-partic le sp ectrum as N increases (as already 
pointed out in Ref.[lO( in the 2D case). This qualitative 
behavior is also valid for the axial symmetric case, where 
two integer frequencies are equal. 

In contrast, in the limit K T > 1 (A < N c ), E(N) is 
controlled by the triaxial term of the shell energy, E T . As 
in the irrational case, the fluctuating part of the energy 
is now characterized by a roughly constant (and small, 
compared to the isotropic case) amplitude, with the gross 
features of the oscillations well approximated by the su- 
perposition of a small number of terms having different 
frequencies. In particular, for some triplet of frequencies, 
in the regime N < N c supershell-like structures may be 
observed, as illustrated for instance in Fig. [3] 

In summary, in the triaxial rational case generically a 
transition at 1Z T ~ 1 (N ~ N c ) from a triaxial regime 
(1Z T 3> 1, N <C N c ) toward an isotropic behavior (1Z T <C 
1, N 3> N c ) is expected as N increases. This behavior 
is illustrated for a particular set of frequencies in Fig. [3j 
where the transition 1Z T ~ 1, estimated from Eqs.(|10p- 
(TT2"]). occurs at A c ~ 131060. The lower left and right 
panels of Fig. [3] focus on the two extreme limits. 

Though we have not studied all the possible cases (ax- 
ial symmetry, rational-irrational, etc), the above exam- 
ples illustrate the variety of behaviors of energy fluctua- 
tions one can find as the frequency ratios are varied. 



or atomic nuclei, behave differently. In those systems the 
shape of the self-consistent mean field potential, deter- 
mined from the minimization of the ground state energy 
with respect to deformations, may strongly depend on 
the number of particles A. 

If the HO potential considered here is viewed as the 
self-consistent mean field potential, then at each A the 
energy is minimized with respect to the three frequencies. 
Conservation of the volume implies that the product of 
the three frequencies remains constant 0], Q1Q2Q3 — 
cte. The constant can be in fact a function of N, to 
mimic for instance the approximately constant nuclear 
density. But this acts as a multiplicative factor which 
can be absorbed into the frequencies. Hence, without 
loss of generality, we normalize this constant to one. 

In the following we consider only axial symmetric de- 
formations (Qi — Q2)- This simplification, though not 
exact, is known to provide a good approximation in most 
nuclear cases. We use the axially symmetric expression of 
the oscillating part of the energy computed analogously 
to the previous sections, with Qi and Q3 two incommen- 
surable real numbers with pairwise irrational ratio, and 
imposing the volume conservation. There is thus only 
one free parameter (Qi), and we obtain 



k=l 



2ir 2 f-f /c 2 sin(7rfc/Q3) 

cos(7r fc/Q?) cos((27rfc(37V) 1 / 3 /Qi)) 

fc=i 



1 °° 



k 2 sin 2 (wk/Ql) 



+ 



^ (-l) fc + 1 cos((2 7 rfc(3A)V 3 g 2 )) 



1 



^ 1 fc=l 



k 2 sin 2 (7rfcQ 3 ) 



Using similar arguments as in previous sections, we keep 
only the first term in each sum of Eq. (fl3| . Besides, as 
we must minimize the energy with respect to Q\, the 
latter is a function of N . The term which comes from 
the first sum has the largest amplitude (this will be con- 
firmed a posteriori, when the A-dependence of Q\ will 
be know) (l5j . So we can remove all other contributions 
from Eq. (|13p . keeping the simple approximation 



2tt z sm.[Tr/Ql) 



(14) 



III. SELF-BOUND SYSTEMS 

Up to now we have considered non-interacting 
fermions in HO potentials that keep their shape (frequen- 
cies) fixed as the number of fermions is varied. Such a 
scheme is a good approximation for dilute Fermi gases in 
optical traps 0, or in quantum dots, where the external 
potential dominates over the mean field interaction en- 
ergy 0- But self-bound systems, like metallic clusters 



For the minimization, we need also to evaluate E(N). 
The approximated (leading-order) TF chemical potential 
Eq. (0) used until now is not sufficiently accurate because 
it leads to a wrong TF ground state: the minimum of 
E doesn't give an isotropic shape. To obtain a better 
description of E we must improve the calculation of fi N . 
This implies solving the full cubic equation given by the 
canonical condition Af(fL N ) = N. We have done that 
using Cardano's formula [16(, and obtained the result 
presented in the third row of Table I. Then, the analytic 



6 




1500 



TV 



FIG. 4: (Color online) Upper panel: Oscillating part of the 
ground state energy as a function of TV computed numeri- 
cally by the minimization of the energy of the free triaxial 
HO (dots). In full line we plot Eq. (|14|l using the phenomeno- 
logical expression (I21|l for s. In dashed line, analytic result 
for Esb valid at mid-shells, Eq. (|20|) . In dotted line, analytic 
result for E S b for the magic numbers (Table II, case B, last 
row). 

Lower panel: Qi = (1 + e) _1,/3 as a function of TV computed 
numerically by the minimization of the energy of the free tri- 
axial HO (dots). In full line we plot Q\ using the phenomeno- 
logical expression (|21|l for e. In dashed line, the mid-shell 
approximation, Eq. (|19[) . 



expression for the energy to be minimized is 

(3TV) 4 / 3 (3TV) 2 / 3 



E(N) 
2(12)2 



24 



-2/3 {2 + 6 2 ) 



(2 + ^ + ^ 1/3 t (2 T (3 r )1/3) ,(15) 



27r 2 sin(7r6 



where 



i/Qi 



K 3 



(16) 



We assume now that the frequency ratio has the per- 
turbative form 9 = 1 + e, where |e| <C 1 (small de- 
formations). Since we know that close to magic numbers 
the system is spherical, we are interested in the behavior 
in-between shells. We thus concentrate our analysis of 
shell effects for values of TV close to the middle of the 
shells. From Eq. p^)) . we find that the middle of the 
shells correspond approximatively to values of TV given 
by (3TV) 1 / 3 = (2i + l)/2 , where i S N. Hence for these 
values of TV, using a Taylor expansion of 6* 1 / 3 and the 
addition formula of the sine function, we get for the fluc- 
tuating part of the energy: 



Esb(N) 



2tt 3 £ 



(17) 



The latter equation is the leading order in e of Eq.fTf 
(the smooth part is of lower order) . The minimization of 



the energy is thus equivalent to the minimization of the 
function (fl~7|) . We cancel the derivative of Eq. (fl~7|) . to get 



2tt(37V) 1 / 3 £ 



tan 



2tt(3TV) 1 / 3 £ 



(18) 



Note that — e is also a solution of Eq. (fl~8|) (prolate-oblate 
symmetry) . In spite of the fact that Eq. (fT5|) has many so- 
lutions, we keep of course the one that gives the smallest 
value of the function (jTTJ) . We have solved numerically 
this equation for e, and find that the solution is well ap- 
proximated by 



e(TV) = ±1.49TV~ 1/3 , 



(19) 



which gives, for the self-bound energy fluctuations at 
mid-shell 



E SB (N) = -0.0152TV 2/3 



(20) 



We have tested numerically these results (cf Fig. 0|. 
They give a very good approximation of E SB (N) in the 
mid-shell region, showing that the axially symmetric 
deformation is a good approximation. The weight of 
the fluctuating term at mid-shell is thus proportional 
to TV 2 / 3 , with a negative slope. This has the same TV- 
dependence as the amplitude of the fluctuating part of an 
isotropic HO (whereas an axial-symmetric fixed poten- 
tial yields instead a TV 1 / 3 dependence), or to the sublead- 
ing, surface term of the smooth part. This enhancement 
shows the fundamental importance of the TV-dependence 
of the frequencies. 

Beyond this simple approximation, we have also found 
a phenomenological expression of e which turns out to 
be a quite good approximation of the deformation for 
arbitrary TV (see bottom part of Fig. 2]) 



e(N) = 



2(3TV)!/3 



arctan 



2tt((3TV) 1/3 - [(37V) 



1/3 



1/2]) 
(21) 

where [a] denotes here the integer part of a. This equa- 
tion shows that the deformation is, up to a global damp- 
ing factor oc TV -1 / 3 , a strictly periodic function of the 
variable x = (3TV) 1 / 3 , of period Ax = 1 (the same is valid 
for the shell correction energy) . It gives an extremely ac- 
curate description of the TV dependence of the frequency, 
Qi = [1 + e(TV)] -1 / 3 , see Fig. [4] Increasing the particle 
number TV from the value of a magic number (correspond- 
ing to Qi = 1), the frequency initially increases, Q\ > 1, 
corresponding to a prolate shape (accordingly, Q 3 < 1). 
At mid-shell, the shape suddenly changes from prolate, 
Qi > 1, to oblate, Qi < 1, keeping its deviation from 
isotropy \Qi — 1| constant. As TV further increases, Q\ 
diminishes, to arrive finally at Q\ — 1 again at the next 
magic number. The cycle starts again, with a decreasing 
overall amplitude. 

Replacing Eq. (f2"T|) in Eq. fT^ij) we obtain a rather good 
approximation of E SB (TV) , see Fig. SJ Using the approx- 



imation e«l, E SB (N) can be written 



(3iV) 2 / 3 sin ( 2tt (3iV) 1/3 (l + e(N)) 1 / 3 



7T 4 arctan 2tt((3V) 1 /3 _ [(3N) 1 / 3 + 1/2]) 



Esb(N) 



(22) 

where e(N) is given by Eq. (|2"Tj) . The results obtained 
with this equation are very similar to those of Eq.fTTJJ, 
shown in Fig. [4j This equation clearly shows the scaling 
TV 2 / 3 of E SB (N). Both the numerator and the denomi- 
nator in Eq. (f22|) are functions of N that oscillate around 
0. The completely different behavior of E SB (N) (a neg- 
ative decreasing function, see Fig. |4|) comes from a deli- 
cate balance between both oscillatory functions. Close to 
magic numbers, the amplitude of Eq. (|22p . not shown in 
the figure, is smaller than the numerical results; further 
corrections need to be included to improve in this region. 

We can check that the amplitude Ai of Eq. (|22[) (first 
term of the first sum in Eq. (|13|) ) is dominant with respect 
to the amplitudes A2 and A3 of the first terms of the two 
remaining sums in Eq. (|13p . respectively. Using Q^ 3 = 
9 = 1 + e, the value of e given by Eq.(|19p. and keeping 
only the first term of a Taylor expansion in e of the sine 
function in the denominators, one can easily check that 
A2/A1 = A3/A1 « 0.07. This justifies the approximation 
used. 

Because the study of E(N) for self-bound system is a 
minimization problem, we are not able to compute an 
explicit and general expression for the average of the 
fluctuating part of the energy, analogous to Eq.® but 
now valid for self-bound systems. Nevertheless, it fol- 
lows from the previous results that the average scales as 
N 2 / 3 . A numerical fit of the average of E SB , denoting 



(E aB (N)) t 



uN 2 ' 3 , 



(23) 



gives a } 



0.0262. 



As an application of the previous results, we consider 
a schematic model of the atomic nucleus, with N non- 
interacting neutrons and N uncharged non-interacting 
protons whose mean-field potential is assumed to have 
an HO shape. In order to mimic the saturation proper- 
ties of nuclear forces and get dimensional quantities we 
multiply E SB by the factor 82A~ 1/3 MeV, where A = 2N 
is the mass number [2j . The previous factor modifies the 
N dependence of the amplitude of E SB , leading to an am- 
plitude now proportional to N 1 / 3 , instead of N 2 / 3 . The 
energy is minimized as previously. 

We compare the results of this model to one of the 
best theoretical models for the nuclear binding energy, 
based on an extension of the liquid drop model Q, see 
Fig. The smooth part used in Ref. [8| being differ- 
ent from ours, it leads to a difference in the offset. In 
spite of that, we see that the qualitative properties of 
the simplified HO model, including the energy scale, are 



correct. In both cases we plot the numerical average of 
E SB . A negative slope of the fit in Ref. [§], consistent 
with — CLsyjN 1 / 3 , is observed, with a proportionality fac- 
tor a sw — 1.93 MeV. It is very similar from what we 
obtain from our simplified model (clho = 1.71 MeV, top 
part of Fig. [5]). The shell amplitude as a function of N 
is also well reproduced. Nevertheless no supershell struc- 
ture occurs in our analysis [l"?] ]. and magic numbers are 
badly estimated because we are not taking into account, 
in our simple model, spin-orbit effects Q. 
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FIG. 5: (Color online) Upper panel: Numerical computation 
of Esb as a function of the neutron number for the HO atomic 
nucleus model (dots). In dashed line, the numerical average 
of E S b and in full line the phenomenological result adapted 
from Eq.JHSJ): {E SB (N)) N = - 1.71 V^MeV 
Lower panel: Ground state shell correction energy as a func- 
tion of the neutron number from (dots). In dashed line, 
the numerical average of E SB and in full line the phenomeno- 
logical fit: {E SB (N)) N = (8.14- 1.93V 1/3 ) MeV. 



IV. CONCLUSION 

For a Fermi gas treated in the mean-field approxi- 
mation, we provided here a quantitative description of 
the ground-state shell correction energy. Two different 
possibilities were considered, fixed potentials and self- 
bound systems. We showed that very different qualita- 
tive behaviors may appear in harmonic traps, depending 
on the number theoretic properties of the frequency ra- 
tios. Based on a semiclassical theory, we provided an 
analytic description of the different behaviors, as well as 
of the average of the shell energy. In spite of its sim- 
plicity, HO potentials show good qualitative agreement 
with more elaborate models. In self-bound systems, one 
of our main results is an accurate analytic description of 



8 



the deformation and of the shell correction energy as a 
function of the particle number. 

A few periodic orbits were used to describe E(N). Al- 
though this is in general a goo d ap proximation, in some 
cases it fails. As discussed in [ll| for jo, the triaxial ir- 
rational case with frequencies very close to each other 
can reproduce, to a good approximation, the isotropic 
case, but only if a large number of terms in the sums 
is included. Gencrically, the convergence properties of 
the sums over periodic orbits is a delicate problem, di- 
rectly related to the number-theoretic properties of the 
frequency ratios. We haven't considered this problem in 
its full generality; it clearly deserves a closer inspection, 
in particular in connection with the minimization prob- 
lem. 

We have shown that, in self-bound systems, the prop- 
erties of the shell correction energy scale as — aJV 2 / 3 (or 
—aN 1 / 3 , if saturation properties are taken into account), 
where a is some positive constant. Although it com- 
pares favorably with realistic nuclear models, it should be 
noted that this behavior is specific of HO potentials. In 
fact, if one considers instead a hard-wall spherical poten- 
tial, supershell structures will appear, and modulations 
of, e.g., the average of E SB (N) are expected, as observed 
in metallic clusters [3]. These, however, are not ob- 
served experimentally in nuclear data, due to the limited 
number of nucleons (see, for instance, 17[). Deviations 



from an HO confining potential toward a more steeper 
one of the Woods-Saxon type are generically expected 
to be produced by interactions. Within a spherical sym- 
metry, these deviations from an harmonic confinement 
were shown to lead to supershell effects [19(. We would 
like to thank Peter Schuck for usefull discussions. 

This work was supported by grants ANR-05-Nano- 
008-02 and ANR-NT05-2-42103, by the IFRAF Insti- 
tute. 



[1] R. B. Wiringa, Steven C. Pieper Phys.Rev.Lett. 89 
(2002) 182501; Steven C. Pieper, K. Varga, R. B. 
Wiringa, Phys. Rev. C 66 (2002) 044310; S. Bargi, G. 
M. Kavoulakis, S. M. Reimann, Phys. Rev. A 73 (2006) 



033613. 

[2] P. Ring and P. Schuck, The Nuclear Many-Body Problem 

(Springer- Verlag, Berlin, 1980). 
[3] A. Bohr and B. R. Mottelson, Nuclear Structure Vol.1 

(Benjamin, Reading, Massachusetts, 1969). 
[4] M. Brack and R. K. Bhaduri, Semiclassical Physics 

(Addison- Wesley, Reading, MA, 1997). 
[5] M. Centelles, P. Leboeuf, A. Monastra, J. Roccia, P. 

Schuck and X. Vinas, Phys. Rev. C 74, 034332 (2006). 
[6] K. M. O'Hara et al, Science 298, 2179 (2002). 
[7] R. Grimm, M. Weidenriiller, Y. B. Ovchinnikov, Ad- 
vances in Atomic, Molecular and Optical Physics 42, 95 

(2000); H. Heiselberg and B. Mottelson, Phys. Rev. Lett. 

88, 190401 (2002). 
[8] P.M. Moller, JR. Nix, W.D. Meyers, W.J. Swiatecki, At. 

Data Nucl. Data Tables 59, 185 (1995). 
[9] S. G. Nilsson, Kgl. Dan. Viden. Selsk. Mat. Fys. Medd. 

29 No.16 (1955) 29 
[10] A. Pandey, O. Bohigas, M. J. Giannoni, J. Phys. A 22, 

4083 (1989). 

[11] M. Brack, S. Jain, Phys. Rev. A 51, 3462 (1995). 
[12] P. Meier, M. Brack, S. C. Creagh, Z. Phys. D 41, 281 
(1997). 

[13] P. Leboeuf and A. G. Monastra, Ann. Phys. 297, 127 
(2002). 

[14] H. Nishioka, K. Hansen, and B. R. Mottelson, Phys. Rev. 
B 42, 9377 (1990). 

[15] We may have chosen a rational frequency ratio between 
Qi and Qi, instead of an irrational one. However, it can 
be shown that the isotropic term that appears for a ratio- 
nal ratio, discussed in section II. C, is not dominant at the 
middle of the shells, even in the limit N — > 00, and that 
the analysis based on rational or irrational frequency ra- 
tios give similar results. This is due to the A r -dependence 
of Qi. 

[16] M. Abramowitz and I. A. Stegun, Handbook of Mathe- 
matical Functions with Formulas, Graphs, and Mathe- 
matical Tables, 9th printing (New York: Dover, 1972). 

[17] P. Leboeuf, Regularity and chaos in the nuclear masses, 
Lectures delivered at the VIII Hispalensis International 
Summer School, Sevilla, Spain, June 2003 (to appear in 
Lecture Notes in Physics, Springer- Verlag, Eds. J. M. 
Arias and M. Lozano); [nucl-th/0406064| 

[18] W. D. Knight, K. Clemenger, W. A. de Heer, W. A. 
Saunders, M. Y. Chou, and M. L. Cohen, Phys. Rev. 
Lett. 52, 2141 (1984). 

[19] YYu, M. Ogren, S. Aberg, S. M. Reimann, M. Brack, 
Phys. Rev. A 72, 051602(R) (2005). 



9 



Q1Q2Q3 



2 _ Qi + gj + gj 
12 



3Q1Q2Q3 



3 gj + Q2 + Qi 

M : M 



S(JV) 



1 n n ah 4 / 3 _i_ Qi + Q2 + Qi /,n n n ah 2 / 3 _l (Qi + Qi + Qi) 2 , n(AT~ 2 / z \ 
QiQ 2 Q 3 V4 ( QlQaQ3 } + 24 (3Q1Q2Q3AO + ^ + 0(N ) 



TABLE I: The smooth part of the single-particle level density for a triplet of frequencies of the 3D harmonic potential as a 
function of the energy e in the first row. Corresponding integrated level density up to an energy fi in the second row. Last row, 
the smooth part of the ground-state energy of the Fermi gas as a function of the number of particles N. 
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A. Triaxial case with irrational frequencies 



B. Isotropic case 



C. Triaxial case with rational frequencies 



-f- 



-l) t+1 cos(27rk/Qi 



Qi sm(TrkQ 2 /Qi) sm(TrkQ 3 /Qi 



2e 2 ^(-l)' t cos(27rfe) 



?i(e) + /5r( E ) 



-y- 

O, ^ si 



-lj'+'cosprfe/Q, 



^) = ^B-i)* Wl+<5a+Qs) 

~ M = J_ V (-l) t+1 cos(2rfe/Q 1 ) 

Oi ^ sin^/ft) sin(jrfcQ 3 /Qi) 
1 A (-l) t,+1 cos(2;rfa/Q 2 ) 
Q 2 ^ sin(7riQ 1 /(3 2 ) S in(rfQ3/(32) 



Q 2 ^ sm{nkQ 1 IQ 2 )sm(iTkQ 3 IQ2 



-y- 

sin 



-l)* +1 cos(2mfce/Q3 



Qs ^ sin(jrfcQi/Q3) sm(7rfeQ 2 /Q 3 ) 



— V - 



(-l) fe+1 co S (27rk/g 3 ) 



sin(7r/;Qi/Q 3 )sin(jr/;(2 2 /Q 3 ) 
=i<-h 



1 A (-l) t+1 sin(27ri:fi/Qi) 
2-k j-> isin(7riQ 2 /Qi)siii(irfeQ 3 /gi) 



-EH) 

7T 



t sin(27rA:^) 



-y- 



(-l) t+1 sin(2^^/Q 2 ) 
2tt ^ fcsin(jr/cQi/Q 2 )sin(jrfcQ 3 /Q 2 ) 

1 y, (-l) t+1 5 m(27rfc M /Q 2 ) 
2jt ^ fesin(jrfcQi/Q 2 ) sin(7rfeQ 3 /Q 2 ) 



Hi{n) + Mt{h) 



_ 1) KQ 1 +fe+Q3)^ sin(2;r ^ ) 
tt/;QiQ 2 Q 3 



MM=y 

" 1 

Qi 

J_ A (-l) t+1 sin(27r^/Q 2 ) 
•2tt 2s 



27r fcsin(7rfcQ 2 /Qi) sin(7r&Q 3 /Qi) 



27r Jzk ksm(irkQi/Q2)^{^kQ 3 /Q 2 ) 
if- (-l) t+1 S in(27r^/Q 3 ) 



£(JV) 



4--' 
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(^l) t+1 CQ5(27rfc(3Q 2 Q 3 A f /Q;) 1 ' 



sin(7ri:Q 2 /Qi) sm(irkQ 3 /Qi) 
(-l) t+1 co 5 (2rf(3QiQ 3 A f /Qi) 1/3 ) 

fc 2 sin(7riQi/Q 2 ) sin(7r/rQ 3 /Q 2 ) 
(-!)"+' cos(2rf(3QiQ 2 jV/Qj) 1/3 ) 
4?r 2 ^ it 2 sin(7rfcQi/Q 3 ) sm{TrkQ 2 /Q 3 ) 
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E 
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^(«) = E 



E,(N) + Et(N) 
(_ 1 ^(qi+q 2 +q 3 )( 3 j V )2/ 3 (3QiQ 2 Q 3 iV) 1/3 
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5 IN) = 9L V (-i) t+1 ^(^(3Q2Q 3 A r /Q;) 1/3 ) 

fc 2 ™(rf<32/Qi)sin( Jr feQ 3 /Qi) 
Qp_ y (-l) t+1 cos(2rf(3QiQ 3 A'/Q|) 1/3 ) 

A (-l)'' +1 co5(2ri(3Q 1 Q 2 JV/Ql) 1 ^ 3 ) 

+ 4jr 2 2_, fc 2 s in(7rfcQi/Q 3 ) S in(vrfc0 2 /0 3 ) 
«ytQ3 



<B(N))a 



& MAGIC 

(JV) 



24 



(3W) ' 



(3N) 2 ' 



24(Q!Q 2 Q 3 )i/3 



for K T <C 1 



24 



(3iVf 



aJV" 273 , forR T >4 
-(3JV) 2/3 



TABLE II: The three columns separate the different cases of the 3D harmonic potential studied in section II. The oscillating 
par t of the single-particle level density for a triplet of frequencies as a function of the energy e is given in the first row (from 
[111]). Corresponding integrated level density up to an energy \i in the second row. Third row, the oscillating part of the 
ground-state energy of the Fermi gas as a function of the number of particles TV. The fourth row gives the analytic expression 
for the average of E. In the last, when available, the analytic curve for the envelope of the magic numbers. 



